subroutine setup_onebody_contributions(it,ip,ij, s_box, gmatrix_configs)
  USE single_particle_orbits
  USE constants
  USE configurations
  IMPLICIT NONE
  TYPE (configuration_descriptor), INTENT(in)  :: gmatrix_configs
  INTEGER :: i,j, ijd, it, ip, ij, ia, ib, ic,id,iabcd,ie,iab, icd, itemp
  REAL(DP), DIMENSION(gmatrix_configs%number_confs, &
       gmatrix_configs%number_confs,n_startenergy_veff), INTENT(INOUT)  :: s_box
  INTEGER :: iph, n_confs, a,b,c,d,na,la,ja,tza,nb,lb,jb,tzb,nc,lc,jc,tzc,nd,ld,jd,tzd, i1, j1, k
  REAL(DP), dimension(n_startenergy_veff) :: onea1,onea2,onea3,onea4, oneb1,oneb2,oneb3,oneb4,s_dir,s_ex
  REAL(DP):: dij, particle_mass
 write(*,*) "hello-li-8"
  DO i1=1,gmatrix_configs%number_confs


     a=gmatrix_configs%config_ab(i1*2-1)
     b=gmatrix_configs%config_ab(i1*2)
     na=all_orbit%nn(a) 
     la=all_orbit%ll(a)
     ja=all_orbit%jj(a)
     tza = all_orbit%itzp(a)
     nb=all_orbit%nn(b)
     lb=all_orbit%ll(b)
     jb=all_orbit%jj(b)
     tzb = all_orbit%itzp(b)
     DO j1=1,gmatrix_configs%number_confs
        c=gmatrix_configs%config_ab(j1+j1-1)
        d=gmatrix_configs%config_ab(j1+j1)
        nc=all_orbit%nn(c) 
        lc=all_orbit%ll(c)
        jc=all_orbit%jj(c)
        tzc = all_orbit%itzp(c)
        nd=all_orbit%nn(d)
        ld=all_orbit%ll(d)
        jd=all_orbit%jj(d)
        tzd = all_orbit%itzp(d)

!        ekin1 = 0.0_dp; ekin2 = 0.0_dp; ekin3 = 0.0_dp ; ekin4 = 0.0_dp
!        upot1 = 0.0_dp; upot2 = 0.0_dp; upot3 = 0.0_dp ; upot4 = 0.0_dp
!        onebd1 = 0.0_dp; onebd2 = 0.0_dp; onebd3 = 0.0_dp ; onebd4 = 0.0_dp   
        

        onea1=0.0_dpc; oneb1=0.0_dpc
        onea2=0.0_dpc; oneb2=0.0_dpc
        onea3=0.0_dpc; oneb3=0.0_dpc
        onea4=0.0_dpc; oneb4=0.0_dpc
        IF ( b == d .AND. la == lc .AND. ja == jc .AND. tza == tzc ) THEN
         call diag_3p1h(a,b,c,d,ij,onea1)
         call diag_4p2h(a,b,c,d,ij,oneb1)
       ENDIF
        IF ( a == c .AND. lb == ld .AND. jb == jd .AND. tzb == tzd ) THEN
         call diag_3p1h(b,a,d,c,ij,onea2)
         call diag_4p2h(b,a,d,c,ij,oneb2)
       ENDIF
        IF ( b == c .AND. la == ld .AND. ja == jd .AND. tza == tzd ) THEN
         call diag_3p1h(a,b,d,c,ij,onea3)
         call diag_4p2h(a,b,d,c,ij,oneb3)
       ENDIF
        IF ( a == d .AND. lb == lc .AND. jb == jc .AND. tzb == tzc ) THEN
         call diag_3p1h(b,a,c,d,ij,onea4)
         call diag_4p2h(b,a,c,d,ij,oneb4)
        ENDIF
!         s_dir=(oneb1+oneb2)/dij(a,b)/dij(c,d)
!         s_ex =(oneb3+oneb4)/dij(a,b)/dij(c,d)
         s_dir=(onea1+onea2+oneb1+oneb2)/dij(a,b)/dij(c,d)
         s_ex =(onea3+onea4+oneb3+oneb4)/dij(a,b)/dij(c,d)
        IF( it /= 0 ) s_box(i1,j1,:) = ( s_dir(:) -iph( (2*ij-jc-jd)/2)*s_ex(:) )
        IF( it == 0 ) s_box(i1,j1,:) = (onea1(:)+onea2(:)+oneb1(:)+oneb2(:) ) 
        write(6,*)a,b,c,d,s_box(i1,j1,n_startenergy_veff/2+1),s_dir(n_startenergy_veff/2+1),s_ex(n_startenergy_veff/2+1)
!        write(*,*)a,b,c,d,s_box(i1,j1,n_startenergy_veff/2+1)
!        stop
    ENDDO
  ENDDO
end subroutine


 SUBROUTINE diag_3p1h(a,b,c,d,jt,onebody_diagram_2)
   USE constants
   USE single_particle_orbits
   IMPLICIT NONE
   INTEGER, INTENT(IN) :: a, c,b,d,jt
   INTEGER ::  j_min, j_max, jtot, h, p1, p2, i, nshell1, nshell2, iosc
   REAL(DP) :: factr
   REAL(DP) ::  val1, val2
   REAL(DP), DIMENSION(n_startenergy_veff) :: w , de 
   REAL(DP), DIMENSION(n_startenergy_veff), INTENT(OUT) :: onebody_diagram_2
   REAL(DP), DIMENSION(n_startenergy_g) :: ans1, ans2
   LOGICAL DENCHECK


   onebody_diagram_2=0.0_dpc
   DO h=1, all_orbit%total_orbits
      IF (all_orbit%orbit_status(h) /= 'hole') CYCLE         
      nshell1=all_orbit%nshell(h)+all_orbit%nshell(c)
      j_min=ABS((all_orbit%jj(a)-all_orbit%jj(h))/2)
      j_max=(all_orbit%jj(a)+all_orbit%jj(h))/2
      DO jtot=j_min,j_max
         factr=(2.*jtot+1.)/(all_orbit%jj(a)+1.)
         DO p1=1, all_orbit%total_orbits
            IF (all_orbit%orbit_status(p1) == 'hole') CYCLE         
            DO p2=1, all_orbit%total_orbits
               IF (all_orbit%orbit_status(p2) == 'hole') CYCLE         
               nshell2=all_orbit%nshell(p2)+all_orbit%nshell(p1)
               iosc=nshell2-nshell1
               IF(dencheck(iosc))CYCLE
              de=all_orbit%evalence(c)+all_orbit%e(h)- &
                   all_orbit%e(p1)-all_orbit%e(p2)+wcn

!             de=all_orbit%e(d)+all_orbit%e(c)+all_orbit%e(h)- &
!                    all_orbit%e(p1)-all_orbit%e(p2)-all_orbit%e(d)+wcn
!               de=starting_energy+all_orbit%e(h)- &
!                    all_orbit%e(p1)-all_orbit%e(p2)-all_orbit%e(d)+wcn
               CALL pphhmtx(a,h,p1,p2,jtot,ans1) 
                IF ( ans1(1) == 0.0_dp) CYCLE
               CALL pphhmtx(p1,p2,c,h,jtot,ans2) 
                 IF (ans2(1) == 0.0_dp) CYCLE
               w=all_orbit%evalence(c)+all_orbit%e(h)+wcn
               DO  i=1, n_startenergy_veff
                  CALL interpolate(w(i),e_start_g,ans1,val1)
                  CALL interpolate(w(i),e_start_g,ans2,val2)
                  onebody_diagram_2(i)=onebody_diagram_2(i)+ &
                       factr*0.5_dp*val1*val2/de(i)
               ENDDO
            ENDDO
         ENDDO
      ENDDO
   ENDDO
 
 END SUBROUTINE
 !
 !     2h1p diagram
 !
    SUBROUTINE diag_4p2h(a,b,c,d,jt,onebody_diagram_3)
      USE constants
      USE single_particle_orbits
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: a, c,b,d,jt
      INTEGER ::  j_min, j_max, jtot, p, i, h1, h2, nshell1, nshell2, iosc
      REAL(DP) :: factr
      REAL(DP):: val1, val2
      REAL(DP), DIMENSION(n_startenergy_veff) :: w1, w2 , de 
      REAL(DP), DIMENSION(n_startenergy_veff), INTENT(OUT) :: onebody_diagram_3
      REAL(DP), DIMENSION(n_startenergy_g) :: ans1, ans2
      LOGICAL DENCHECK
    

      onebody_diagram_3=0.0_dp
      DO p=1, all_orbit%total_orbits
         IF (all_orbit%orbit_status(p) == 'hole') CYCLE         
         j_min=ABS((all_orbit%jj(a)-all_orbit%jj(p))/2)
         j_max=(all_orbit%jj(a)+all_orbit%jj(p))/2
         DO jtot=j_min,j_max
            factr=(2.*jtot+1.)/(all_orbit%jj(a)+1.)
            DO h1=1,all_orbit%total_orbits
               IF (all_orbit%orbit_status(h1) /= 'hole') CYCLE         
               DO h2=1,all_orbit%total_orbits
                  IF (all_orbit%orbit_status(h2) /= 'hole') CYCLE         
                  nshell1=all_orbit%nshell(h1)+all_orbit%nshell(h2)
                  nshell2=all_orbit%nshell(p)+all_orbit%nshell(a)
                  iosc=nshell2-nshell1
                  IF(dencheck(iosc)) CYCLE
              de=all_orbit%e(h1)+all_orbit%e(h2)- &
                   all_orbit%e(p)-all_orbit%evalence(a)+wcn

!                  de=all_orbit%e(d)+all_orbit%e(c)+all_orbit%e(h1)+all_orbit%e(h2)- &
!                       all_orbit%e(p)-all_orbit%e(a)-all_orbit%e(d)-all_orbit%e(c)+wcn
                  CALL pphhmtx(h1,h2,c,p,jtot,ans1) ; IF ( ans1(1) == 0.0_dp) CYCLE
                  CALL pphhmtx(a,p,h1,h2,jtot,ans2) ; IF (ans2(1) == 0.0_dp) CYCLE
                  w1=all_orbit%e(h1)+all_orbit%e(h2)+ &
                       all_orbit%evalence(c)-all_orbit%evalence(a)+wcn
                  w2=all_orbit%e(h1)+all_orbit%e(h2)+wcn
                  DO i=1, n_startenergy_veff
                     CALL  interpolate(w1(i),e_start_g,ans1,val1)
                     CALL  interpolate(w2(i),e_start_g,ans2,val2)
                     onebody_diagram_3(i)=onebody_diagram_3(i)- & 
                          factr*0.5_dp*val1*val2/de(i)
                  ENDDO
               ENDDO
            ENDDO
         ENDDO
      ENDDO
    
    END SUBROUTINE 

